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ABSTRACT 

We perform hydrodynamical simulations of major galaxy mergers using new methods 
for calculating the growth of massive black holes (BH) in galactic nuclei and their im- 
pact on the surrounding galaxy. We model BH growth by including a subgrid model for 
accretion produced by angular momentum transport on unresolved scales. The impact 
of the BHs radiation on surrounding gas is approximated by depositing momentum 
into the ambient gas, which produces an outward force away from the BH. We argue 
that these phenomenological models for BH growth and feedback better approximate 
the interaction between the BH and dense gas in galaxies than previous models. We 
show that this physics leads to self-regulated black hole growth: during the peak of 
activity, the accretion rate onto the BH is largely determined by the physics of BH 
feedback, not the subgrid accretion model. The BH significantly modifies the gas dy- 
namics in the galactic nucleus (< 300 pc), but does not generate large-scale galactic 
outflows. Integrated over an entire galaxy merger, BH feedback has little effect on the 
total number of stars formed, but is crucial for setting the BHs mass. 



1 INTRODUCTION 



Modern theories of galaxy formation hold that strong feed- 
back processes regulate star formation in galaxies across 
a wide range of masses. For more massive galaxies, stel- 
lar feedback processes appear to become less efficient and 
feedback from a central massive black hole (BH) begins to 
dominate. Feedback from an active galactic nucleus (AGN) 
has been invoked to account for many observational results 
in galaxy formation, including the Mbh — and Mbh — o 
relations and the suppression of star formation in elliptical 
galaxies (ISilk fc Re e^'lQQS^: 'King"20o5; ' Murrav et all (20051: 



[Di Mat teo et al 
|2007|)- 



2005; Springcl ct al. 20053; iHopkins et all 



Many recent studies developing numerical models for 
the effects of BHs on galactic scales have used broadly simi- 
lar implementations of the uncertain physics o f AGN fueling 
and feedback (e.g., .S pringcl ct al. 2005b; Ka wata fc GibsonI 
I2OO5I : [Johansson et al.ll2009l '). It is, e.g., often assumed that 
a BH of mass Mbh accretes mass from the surrounding in- 
terst ellar medium (ISM) at a rate proportional to the Bondi 
rate l|Bondilll952l ). 



(1) 



where p is the surrounding gas density, Ca is the sound speed 
and / ~ 10 — 100 takes into account the fact that the 
sphere of influence of t he black hole is not always resolved 
ijBooth fc Schavell2009l ). Moreover, these same calculations 
assume that the black hole's impact on its host galaxy can 
be approximated by depositing thermal energy released by 



accretion back into the surrounding gas. There is, however, 
little detailed motivation or justification for either of these 
assumptions. 

Eq. ([T]) is only applicable when the gas fueling the cen- 
tral BH has very little angular momentum. Otherwise, the 
transport of angular mom entum regulates the a ccretion rate 
onto the central BH (e.g.. lShlosman et al.lll990h . It is gener- 
ally believed that gas-rich disk galaxies are the progenitors 
of today's > L* ellipticals and, in particular, that mergers of 
gas-rich galaxies lead to luminous starbursts and the growth 
of central massive BHs (|Sanders et alj|l988l : [ifopk ins et al.l 
l2005h . Most of the gas in disk galaxies, merging galaxies, lu- 
minous starbursts (jPownes fc Solo mon 1998; Tacconi et "al] 
,2006 ) ■ and nearby luminous AGN (Ho ct al. 2008) appears 
to reside in a rotationally supported disk. There is therefore 
no strong reason to believe that the spherically-symmetric 
Bondi accretion rate is a reasonable estimate of the accretion 
rate onto a BH in gas-rich disk galaxies. 

The energy generated by a central AGN can couple to 
its surroundings in a variety of ways, all of which may have 
a significant dynamical infiuence on gas in the host galaxy 
and in the surrounding intergalactic medium. For example, 
relativistic jets inject energy into intracluster plasma and 
may be the key mechanism for suppressing cooling fiows in 
galaxy clusters dMcNamara fc Nulscn 2007), although the 
details of how the energy in the jet couples to the surround- 
ing plasma in a volume-fillin g way are not fully understood 
IVernaleo fc Revnoldj|2006D . On galactic scales, winds from 
an accretion disk around the BH may sweep-up and drive gas 
out of galaxies (e.g.. iKing,,2003. '). And the AGN's radiation 
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can strongly impact the surrounding gas, bot h via Comp- 
ton coohng/heating (e.g., ISazonov et al.l l2005l) and via the 
momen tum imparted as UV radiation is abso rbed by dust 
grains (|Chang et all 1 19871 : ISanders et all ligSSl ) . 

The precise physical mechanism (s) responsible for AGN 
feedback are not fully understood, particularly on galactic 
scales. For this reason, it is useful to distinguish between 
two classes of models: energy and momentum injection. We 
believe that momentum injection, not energy injection, is 
likely the dominant form of feedback for the majority of the 
gas in a galaxy. In most circumstances, jets take the path 
of least resistance and travel relatively unimpeded out of 
a galaxy. Furthermore, while radiation from an AGN can, 
in principle, Compton heat the surrounding gas enough to 
unbind it, it can only do so for very low density gas. For 
example, for a BH radiating at ~ 10^® erg s~^ with a typical 
quasar spectrum, only gas with n < 1 cm~'^ can be heated 
to the Compton temperature within ~ 100 pc. However, the 
mean gas densities in the central ^ 100 pc of ultraluminou s 
infrared galaxies are ~ 10 cm~ (jPownes fc SolomonllT99i ). 
At these densities, the cooling time of the gas is very short 
and the gas is unable to retain any injected energy. Thus if 
the radiation from a BH strongly modifies the dynamics of 
the gas in its immediate vicinity, it must be via the force 
exerted when the radiation is absorbed. 

Given the uncertainties in the physics of BH accretion 
and feedback, it is important to explore a range of models for 
the impact of BHs on galaxy formation. Towards this end, 
we have carried out numerical simulations of major galaxy 
mergers, qualitatively taking into account the physics of ac- 
cretion induced by angular momentum transport and AGN 
feedback by momentum injection (radiation pressure). Our 
accretion and feedback prescriptions both differ from those 
used in previous numerical simulations of BH growth and 
feedback. The results in this Letter are taken from a larger 
set of calculations (DeBuhr et. al. in prep.) and represent 
general features of all the simulations we have carried out. 

2 METHODS 

W e use a non-pub lic update of the TreeSPH code GADGET- 
2 (|Springe]||2005l l to perform simulations of galaxy mergers 
with feedback from both star formation and central super- 
massive BHs. The code, provided by V. Springel, includes 
the ef fective star formation model of ISpringel fc Hernguistl 
We describe below the additional modifications that 
we have implemented to model BH growth and feedback. 

The multipha se equation of state of 

ISpringel &: Hernguistl (|2003l ) overpredicts the "sound 
speed" as compared to observations of the random veloci- 
ties in galaxies (in atomic or molecular gas). For example, 
the parameter choices Tsn = 4 x 10*K, = 4000, 
tl = 8.4Gyr, a nd qsos = 0.5, whi ch have been used in 
previous works (|Springel et al.|[2005bl ). predicts dv ~ 30 km 
s~^ at n ~ 1 cm~^ and dv ~ 110 km s~^ for n ~ lO^cm"^. 
These are too large by a factor of ^ 2 — 3 compared with 
the observed values (|Downes fc Solomon|[i998 '). To account 
for this difference, we reduce the pressure everywhere by a 
constant factor of 10. As in ISpringel fc Hernguistl l|2003h . 
we assume that p t, ~ p^'^ for consistency with observations 
jKennicuttll 19881 ). The normalization of the star formation 
prescription is chosen so that a Milky Way-like galaxy has 



a total star formation rate of about one solar mass per 
year; for galaxies with different surface densities, the result 
is also consistent with Kcnnicutt (1988). By reducing the 
pressure at fixed p by a factor of 10, the gas is more dense 
in hydrostatic equilibrium. This would increase the star 
formation rate relative to the observed value. To correct 
for this, we modify the equation of state parameters to: 
tl = 13.86 Gyr, /3 = 0.1, Ao = 6600, Tsn = 6.6 x 10* K, 
Tc = 1000 K and qeos = 0.5. 

The simulations described in this work are all mergers 
of equal mass galaxies. Each model galaxy consists of a dark 
matter halo, a rotationally supported disk of gas and stars, 
a stellar bulge, and a central BH. The galax y parameters 
are similar to those in ISpringel et al] (|2005bl ): each galaxy 
has a total mass of 1.36 x IO^^Mq; the mass of the disk 
is 4.1% of the total, i.e., 5.57 x 10^ Mq, where 10% of the 
disk mass is assigned to gas and 90% to stars; the bulge has 
a mass of 1.86 x lO^'' Mq, i.e., 1/3 of the total disk mass. 
Each galajcy is made of 8.0 x 10"" simulation particles and 
the gravitational force softening is e = 47 pc. The halo and 
the bulge have Hernquist ( 1 99Q) density profiles, where the 
virial and half-mass radii of the halo are 229 kpc and 102 
kpc, respectively (the concentration is 9.0), and the effec- 
tive radius of the bulge is 1.27 kpc. The gaseous and stellar 
disks have exponential profiles with scale lengths Rd — 3.51 
kpc; the scale-height of the stellar disk is 0.7 kpc, while the 
scale-height of the gaseous disk is determined by hydrostatic 
equilibrium. The massive BH in each galaxy is modeled us- 
ing a specially marked coUisionless tracer pa rticle. 

Th e initial conditions are generated as in lSpringel et al.l 
(j2005bl ). except for the decrease in gas pressure described 
above. The galaxies are placed on a prograde parabolic orbit. 
The individual spins of the two galaxies are randomly chosen 
to have a relative angle of about 41 degrees. The galaxies 
begin at a distance of 142 kpc and the orbit has a pericenter 
of 14.2 kpc. 

We estimate the accretion rate onto the BH from the 
surrounding gas, due to viscous transport of angular mo- 
mentum, using 

M„i. = SttqE^ , (2) 

where E is the mean surface density of the gas in the disk, Q 
is its angular rotational frequency, and the free parameter a 
is the dimensionless viscosity; M^is is also capped at the Ed- 
dington rate. We compute E and the sound speed Cs by tak- 
ing an average of the properties of the individual SPH par- 
ticles in a spherical region with a radius Race = 4e = 188 pc 
centred on the BH. Using velocity information directly from 
the simulation particles themselves to compute SI proved to 
be too noisy; we thus determined Q, using the total mass, 
Mt, inside Race through Q,'^ = GMt/RLc- 

Note that in our accretion prescription, Mvis = if 
there is no gas within 4e of the BH. This feature of our 
model accounts for the fact that our simulations capture 
the angular momentum transport produced by gravitational 
torques on large scales (> 4e = Race); we assume, as is 
physically reasonable but by no means proven, that these 
must be sufficient to bring gas close to the BH (within Race) 
for significant BH accretion to proceed. 

Eq. ^ is rem i niscen t of an alpha prescription of 
IShakura fcSunvaevI (|l973l ). but in this formulation a pa- 
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rameterises both the efficiency of angular momentum trans- 
port on scales smaller than our gravitational force softening 
(by, e.g., gravitational torques) and the uncertainty due to 
the fraction of the inflowing mass that turns into stars vs. ac- 
creting onto the central BH. The physical processes responsi- 
ble for transportin g gas from ~ kp c to ~ 0.1 pc are still not 
fully understood ( GoodmarJ 120031 ) , but n on-axisymmetric 
gravi tational torques are likely responsible l|Shlosman et al.l 
Il989l ). Detailed calculations of the structure of AGN disks 
from ~ 0.01 — 100 pc, based on transport by spiral waves, 
show that Eq. ((Sjl, evaluated at radii ~ 30 — 100 pc, can 
provide a reasona ble estimate of the accr etion rate onto the 
BH in some cases jThompson et al.ll2005l ). Although Eq. ^ 
is only a crude approximation to the true accretion rate 
onto the BH, it captures the qualitative physics of accre- 
tion induced by angular momentum transport, and is thus, 
we believe, a more suitable "subgrid" model than Eq. ([T]). 
Our fiducial choice for a is a = 0.05, but we also present 
results for 0.15 in Sec. 3. In the future, we intend to better 
calibrate our model of angular momentum transport using 
simulations that focus on t he central ~ 100 pc of galaxies 
iHopkins fc QuataertllioogI ). 

To model the feedback onto the gas surrounding the 
BH, we have implemented a simple phenomenological model 
in which the AGN's luminosity L is coupled back into the 
surrounding gas by depositing momentum radially outward 
from the BH. Our goal is to account for the radiation pres- 
sure produced by the absorption and scattering of the AGN's 
radiation by dust in the ISM. To accurately do so would re- 
quire a radiation transport calculation, which is beyond the 
scope of the current paper. Instead, we model the impact 
of this radiation on the surrounding galaxy by depositing a 
total momentum (per unit time) of 



tL/c, where L = min(r;M„isC , Lfidd) 



(3) 



radially away from the BH into every SPH particle within a 
distance of Race from the BH; each particle receives the same 
acceleration. Note that the number of particles that receive 
this extra force, A'', will change with time as particles enter 
and leave the central region of radius Race- We assume a 
radiative efficiency of 77 — 0.1. Equation ((3| models the ab- 
sorption by the dust of the UV radiation from the AGN (one 
L/c), and, more importantly, the subsequent diffusion of the 
far IR photons (r L/c) . In this way, the value of r determines 
the total momentum deposited and corresponds to the total 
far IR optical depth in the nuclear region; we choose r = 10 
in these calculations. This value of r is consistent with the 
fact that even the far-infrared radiation produced by dust is 
optically thick at radii ~ 100 pc during galaxy mergers (e.g., 
[Thompson et al.ll2005l ). In particular, r ~ 10 is motivated 



by the surface density of E ~ 3 g cm in the inner 



100 



pc near the peak of accretion, and a FIR opacity of ~ 3 cm 
g~^. The exact value of r does not significantly affect our 
conclusions, but it does normalize the values of Mvis and 
Mbh (see Eq. g] below) . 

The strength of the feedback on an individual particle 
depends not only on the luminosity, but on the number of 
particles, A'', to which the force is being applied in a given 
timestep. Our results do not depend strongly on N; this is 
because the momentum is quickly shared with the rest of the 
gas particles via pressure forces. We carried out a number 
of test problems on the evolution of gaseous shells with the 



additional force tL/c; these explicitly show no dependence 
on N (DeBuhr et. al. in prep). 

Computing the accretion rate onto and the feedback 
from the BH in the simulations is prone to noise induced 
by the stochastic motion of the BH particle. To avoid this 
"Brownian" motion, we choose a mass for the BH tracer 
particle of 2.8 x IO^A/q, which is roughly a factor of 100 
higher than the other particle masses in the simulation. Note 
that this mass is an artificial dynamical mass for simulation 
purposes; in addition to this, we integrate M{t) to determine 
the "true" Mbh- Once the two BH tracers have a separation 
of ARacc or smaller, we consider that they would coalesce to 
form a single BH. Once the tracers merge, the two values 
of Mbh are summed and one of the BH particles is moved 
to the center of mass of the two tracers, and the other is 
removed from the region. 

3 RESULTS 

The top panel of Figure [T] shows the viscous accretion rate 
onto the BH for the fiducial run with a — 0.05 (black) and 
for a run with q = 0.15 (blue). For comparison, the Edding- 
ton rate MEdd = LEdd/T? is shown in grey, using the BH 
mass as a function of time from our fiducial simulation. The 
accretion rate is relatively constant at early times but then 
peaks during the first close passage of the two galaxies at 
^ 0.75 Gyr and then even more strongly as the two galaxies 
complete their merger at ~ 1.6 Gyr; note that M <^ MEdd 
at both early and late times but reaches ~ MEdd for ~ 100 
Myrs near both first and final passage. 

One of the interesting results in Figure [1] is that dif- 
ferences in a do not significantly change the accretion rate 
onto the BH, particularly near the peaks of activity. This 
is contrary to what one might expect from the fact that 
Mvis oc a (eq. [2]). The origin of the weak dependence of 
Af„is on a is that when the supply of mass is large, feedback 
from accretion onto the BH regulates the rate at which the 
BH accretes. Previous work has shown that there is a crit- 
ical luminosity Lc at which the outward radiation pressure 
force due to the central AGN just balances the inward force 
of gravity. For a simple spherically symmetric problem, this 
is given by tLc/c = Afga'^/G, where /g is the gas fraction in 
the nuclear region and = GMt /2R acc, with Mt the total 
mass within Race (iMurrav et aLlfioOSl ). This in turn implies 
a critical accretion rate Mc 



Mc 



rricG 



(4) 



For L > Lc, the radiation force on gas in the nuclear re- 
gion exceeds the inward force of gravity, and thus gas in 
the vicinity of the BH will be pushed out of the nuclear 
region. In our model, the accretion rate is determined by 
the gas properties withm Race — 188 pc; thus if the gas is 
largely pushed out of the nuclear region, the accretion rate 
onto the BH decreases. When L < Lc, gravitational torques 
can drive gas into the nuclear region towards the BH, thus 
increasing M„is. This suggests that the accretion rate may 
self-adjust such that M^jis ~ Mc- To quantify this, the top 
panel of Figure [1] shows Ale computed within 2Racc of the 
BH for the fiducial calculation (red). The accretion rate is 
indeed ~ Mc near the peaks of activity. This highlights that 
although Mvis <C Mc is certainly possible if there is insuf- 
ficient gas in the nuclear regions (e.g., after the merger). 
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Figure 1. Top panel: Viscous accretion rate Myis onto tiie black 
hole (eq. [2]) for our fiducial galaxy merger simulation {a = 0.05; 
black), and for a run with three times the viscosity {a = 0.15; 
blue). Also shown is the critical mass accretion rate Mc (eq. [4]) 
at which radiation pressure can push gas out of the nuclear re- 
gion (red), and the Eddington rate M^dd (grey). Note that the 
accretion rate adjusts to Mvig ~ Mc during the peaks of activity, 
independent of a. Bottom panel: Viscous accretion rate for the 
fiducial simulation (black) as compared to a simulation without 
feedback (orange; a = 0.15); also shown is the Bondi accretion 
rate using the BH mass from the fiducial simulation. Unlike Al^ig, 
the Bondi rate is ~ M^dd a-t nearly all times. 



feedback limits the maximum rate at which the BH can ac- 
crete to be ~ Mc- One point that we return to below is that 
although feedback does have a strong effect on the gas dy- 
namics in the galactic nuclei, it is not strong enough to blow 
large amounts of gas out of the galaxy as a whole. 

The bottom panel of Figure \T\ compares the accretion 
rate for the fiducial run (black) and a similar run with no 
BH feedback (orange; a — 0.15). The peak accretion rate is 
a factor of ~ 10 higher in the case without feedback, and the 
duration of activity is significantly longer; moreover, M^is 
in the absence of feedback is oc a and so can be scaled up 
or down by arbitrary amounts by varying a, unlike in the 
presence of feedback (top panel). Also shown in the lower 
panel is the Bondi accretion rate (grey) with the BH mass 
set by that in the fiducial run (which uses M^is). For nearly 
all of the simulation, MBondi — M sdd, and unlike the simu - 
lations of lSpringel et all l|2005bl 'l or [Johansson et all l|2009f ). 
the Bondi rate does not decrease after the final merger, be- 
cause the ambient gas remains cool and dense. 

Figure [2] shows the column density of gas in the vicinity 
of the BH for the case without feedback (right) and the case 
with feedback (left) both at the same time during the peak 
of accretion at t = 1.74Gyr. The images are 28.5kpc on a 
side in the top row, and 4.28kpc on a side in the bottom 



Figure 2. Column densities of gas centered on the BH for the run 
with feedback (left) and without feedback (right), at t = 1.74 Gyr 
during the peak of activity (Fig. [TJ. The colour indicates the col- 
umn density with brighter colour indicating larger densities. The 
large scale images are projected onto the orbital plane, whereas 
the small scale images are viewed along the orbital plane as the 
nuclear disk has a large inclination. The white box shows the rela- 
tive scale of the images. Note that the large-scale gas distribution 
is quite similar in the two cases; the BHs only visible effect is in 
the nuclear region (< 300 pc), where some of the gas is cleared 
out in the case with fe edback. These images were generated using 
SPLASH l |Pricj|2007ll . 



row. In the simulation with feedback, one can see explic- 
itly that the gas has been evacuated from the region near 
the BH (within ~ Race), as argued above. These images 
also demonstrate that the feedback from the BH does not 
produce a large-scale blow-out of matter from the galactic 
nucleus. More quantitatively, at the end of the simulation, 
the runs with and without feedback have the same mass of 
gas outside iRd to within 10%, and the mass that is at large 
radii is due to the merger dynamics (e.g., tidal tails) rather 
than the BH driving a powerful outflow. This is qualitatively 
consistent with observational evidence for large reservoirs 
of atomi c and molecular gas in nearby luminous AGN and 
quasars (IScoville et al.ll2003l ). which have relatively normal 
kinematics ||Ho et al. 20081 ) . By contrast, previous simula- 
tions using the Bondi accretion rate and thermal feedback 
find that the BH unbinds the remaining gas in the galaxy 
near the end of the merger (e.g.,^i Mattco ot al. 2005) and 
that this c an be important for sh utting off star formation in 
ellipticals (|Springel et al.ll2005ah . 

In Figure [3] we show the BH mass Mbh and the in- 
tegrated mass of stars formed during the merger (M,) for 
the fiducial simulation (black), the run without feedback (or- 
ange, dash-dot), and the run with a — 0.15 (blue, dash). The 
addition of feedback changes the total star formation during 
the merger by less than 1%. In addition, the final BH masses 
for the two runs with feedback differ by about 30%. This is 
a further consequence of the self regulated accretion during 
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Figure 3. The black hole mass Mgu (lower three curves) as a 
function of time during the merger, and the mass of new stars 
formed during the merger M, , for three simulations: runs includ- 
ing BH feedback with a = 0.05 (black) and a = 0.15 (blue, 
dashed), and the run without BH feedback (orange, dash-dot). 
Note that BH feedback has little effect on the total stellar mass 
formed, while it reduces the final BH mass by a factor of ~ 10 
compared to the run without feedback. The BH mass is, however, 
nearly independent of a because the accretion self-regulates as 
shown in Figure [l] 



most stages of the merger. There is, of course, freedom in 
choosing the initial mass of the BH in our calculations, but 
so long as this is sufficiently small, it does not significantly 
change the final mass of the BH. For the simulations without 
feedback, the final BH mass is larger than in the presence 
of feedback by a factor of ~ 10, as would be expected from 
Figure [T] In addition, the mass of the BH scales oc a in this 
case. Although we have not made a quantitative compari- 
son, the small dispersion in BH mass for different subgrid 
accretion models in the presence of feedback appears consis- 
tent with the small scatter in the Mbh — o relation (while 
models without feedback would produce a larger dispersion 
in Mbh)- 

4 DISCUSSION 

The model presented here is a necessarily simplified treat- 
ment of the physics occurring in the nuclear regions of galax- 
ies. In particular, the choices of a in the accretion model 
and of Race and r in the model of momentum deposition 
are somewhat uncertain. Changing these values affects some 
of details of the gas dynamics. For instance, the accretion 
history of the BH change modestly as we vary Race and 
r. However, it is encouraging that many of the global re- 
sults of the simulations are insensitive to these choices. The 
peak luminosity occurs at the same time and always reaches 
the Eddington limit. The total stellar mass formed during 
a merger is essentially independent of these parameters and 
the final BH mass is relatively insensitive to both Race and 
Q. Perhaps most interestingly, the peak accretion rate is rel- 
atively independent of the subgrid accretion model (q) and 
is instead set by the structure of the host galaxy and the 
feedback physics, reaching the critical rate ~ Mc at which 
radiation pressure balances gravity in the nuclear regions of 
the galaxy (eq. Q. 

A clear next step in this modeling effort is to perform 
radiative transfer calculations simultaneously with the SPH 



calculation in order to more reliably determine the radia- 
tion pressure force. This would not only eliminate the need 
to specify the parameters Race and t by hand, but would 
also provide information about the AGN spectrum as a func- 
tion of time. Detailed comparisons between these results and 
observations should allow quantitative tests of the impor- 
tance of AGN feedback by momentum deposition during BH 
growth and galaxy formation. 
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